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ABSTRACT 

We present an analytic model for the fully nonlinear power spectrum P and bispec- 
trum Q of the cosmological mass density field. The model is based on physical properties 
of dark matter halos, with the three main model inputs being analytic halo density pro- 
files, halo mass functions, and halo-halo spatial correlations, each of which has been well 
studied in the literature. We demonstrate that this new model can reproduce the power 
spectrum and bispectrum computed from cosmological simulations of both an n = —2 
scale-free model and a low-density cold dark matter model. To enhance the dynamic 
range of these large simulations, we use the synthetic halo replacement technique of Ma 
& Fry (2000), where the original halos with numerically softened cores are replaced by 
synthetic halos of realistic density profiles. At high wavenumbers, our model predicts 
a slope for the nonlinear power spectrum different from the often-used fitting formulas 
in the literature based on the stable clustering assumption. Our model also predicts a 
three-point amplitude Q that is scale dependent, in contrast to the popular hierarchical 
clustering assumption. This model provides a rapid way to compute the mass power 
spectrum and bispectrum over all length scales where the input halo properties are 
valid. It also provides a physical interpretation of the clustering properties of matter in 
the universe. 

Subject headings: cosmology : theory - dark matter - large-scale structure of universe 



1. Introduction 



Two conceptual pictures of galaxy clustering have been examined in the literature, the con- 
tinuous hierarchical clustering model and the power-law cluster model (Peebles 1980, §61). In the 
hierarchical clustering model, which has emerged as the accepted model over the past two decades, 
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galaxy clustering is characterized by power-law correlation functions: the TV-point correlation func- 
tion £at scales with configuration size as £n oc r _7JV oc ^ N *\ where 7^ = (N — 1)7 and the 
two-point correlation function goes as £2 = £ °c r_7 - The hierarchical model is motivated by the 
observed power-law behavior 7 « 1.8 of galaxy correlations (Groth & Peebles 1977; Fry <fe Peebles 
1978), with a theoretical basis in a self-similar, scale-invariant solution to the equations of motion 
(Davis & Peebles 1977). 

The alternative power-law cluster model has an even longer history (Neyman & Scott 1952; 
Peebles 1974, 1980; McClelland & Silk 1977; Scherrer & Bertschinger 1991; Sheth & Jain 1997; 
Valageas 1998; Yano & Gouda 1999). In this model, galaxies are placed in spherical clumps that are 
assumed to follow a power-law density profile p(r) oc r~ e , with the centers of the clumps distributed 
randomly. The resulting two-point correlation function is also a power law with a logarithmic slope 
7 = 2e — 3. While it is possible to reproduce the observed two-point function by an appropriate 
choice of the power index e = (3 + j)/2 « 2.4, Peebles and Groth (1975) pointed out that this 
model produces a three-point function that is too steep to be consistent with observations in the 
Zwicky and Lick catalogs. 

In an earlier paper (Ma & Fry 2000a), we have shown that in the nonlinear regime, the 
three-point correlation function £ = £3 of the cosmological mass density field does not exactly 
follow the prediction £ oc £ 2 of the hierarchical clustering model. These conclusions are drawn 
from study of high resolution numerical simulations of a cold dark matter (CDM) model with 
cosmological constant and of a model with scale-free initial conditions P(k) ~ k n with n = —2. 
In experiments replacing simulation dark matter halos with power-law density profiles, p(r) ~ r~ e , 
we have demonstrated that the behavior of the correlation functions in the nonlinear regime are 
determined by the halo profiles, but that it is not possible to match both the two- and three-point 
correlations with a single slope e. These results differ from the predictions of both of these two 
conceptual models. 

In this paper, we expand our previous study of the nonlinear two- and three-point correlation 
functions by investigating a new prescription that takes into account the non-power-law profiles 
of halos, the distribution of halo masses, and the spatial correlations of halo centers. Each of 
these ingredients has been well studied in the literature. We find that this halo model provides a 
good description of the two- and three-point correlation functions in both the n = —2 and CDM 
simulations over the entire range of scales from the weak clustering, perturbative regime on large 
length scales, to the strongly nonlinear regime on small length scales. Our result is approximately 
hierarchical over an intermediate range of scales, thus uniting the two pictures. An independent 
recent study by Seljak (2000), which appeared during completion of this work, has also examined 
the two-point power spectrum in a similar construction and has found that this type of approach 
can reproduce the power spectrum in the CDM model. The analytic model proposed here can be 
used to compute the two- and three-point correlation functions and their Fourier transforms, the 
power spectrum and bispectrum, over any range of scale where the input halo properties are valid. 
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In a subsequent paper (Ma & Fry 2000c), we study the predictions of this analytic halo model 
for the asymptotic nonlinear behavior of the iV-point correlation functions and the pairwise veloc- 
ities and examine the conditions required for stable clustering. 

The outline of this paper is as follows. In §2 we describe the three input ingredients of the 
model: halo density profiles, halo mass functions, and halo-halo correlations. In §3 we assemble 
these ingredients and construct analytic expressions for the two-point correlation function £(r) and 
the power spectrum P(k). In §4 we do the same for the three-point correlation function C( r i, r 2, ^3) 
and its Fourier transform, the bispectrum B(k\,k2, £3). In §5 we test the validity of this new model 
by comparing its predictions with results from numerical simulations of an n = — 2 scale free model 
and a low-density CDM model with a cosmological constant (ACDM). We also present results of 
the synthetic halo replacement technique used to enhance the numerical resolution. In §6 we discuss 
further the physical meanings and implications of the model. In particular, we elaborate on two 
important implications of this model: deviations from the common assumptions of stable clustering 
and hierarchical clustering. Section 7 is a summary. 



2. Model Ingredients 

2.1. Halo Mass Density Profile 

It has been suggested recently that the mass density profiles of cold dark matter halos have 
a roughly universal shape, generally independent of cosmological parameters (Navarro, Frenk, & 
White 1996, 1997) 

^- = Su(r/R s ), (1) 

where 5 is a dimensionless density amplitude, R s is a characteristic radius, and p is the mean 
background density. We consider two functional forms for the density profiles 

1 



ui(x) = 



xP(l + x) 3 'P 



Unix) = — 7- 5 — r . (2) 

Both forms have asymptotic behaviors x~ p at small x and x~ 3 at large x, but they differ in the 
transition region. The first form uj(x) with p = 1 is found to provide a good fit to simulation halos 
by Navarro et al. (1996, 1997), whereas the second form uji(x) with a steeper inner slope p = 3/2 
is favored by Moore et al. (1998, 1999). Some independent simulations have produced halos that 
are well fit by the shallower p = 1 inner slope (e.g., Hernquist 1990; Dubinski &; Carlberg 1991; 
Huss, Jain, & Steinmetz 1999), and others the steeper p > 1 slope (e.g., Fukushige and Makino 
1997). Jing & Suto (2000) have recently reported a mass-dependent inner slope, with p « 1.5 for 
galactic-mass halos and p « 1 for cluster-mass halos. Many of these authors find that the outer 
profile scales as r -3 , but steeper outer profiles have also been suggested (Hernquist 1990; Dubinski 
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& Carlberg 1991). Given these uncertainties, we will consider in this paper both types of profiles 
in equation (2). 

The parameters R s and 5 in equation (1) are generally functions of the halo mass M. A 
concentration parameter, 

R200 , , 

can be used to quantify the central density of a halo (Navarro et al. 1997), where -R200 is the 
radius within which the average density is 200 times the mean density of the universe. Using 
M = 8007rpi?20o/3> we can relate R s and 5 to M and c, where the scale radius R s is 

1 / 3M \ 1/3 1.63 x 10~ 5 / M \ 1/3 , 

R ' = ~c (soo^j — ig^Hs^) ft Mpc - (4) 

and the density amplitude 5 is 

= 200 c 3 

1 ~ 3[ln(l + c) - c/(l + c)} ' P ~ 1 ' 
, 100 c 3 3 

5/7 = ln(l + C 3/ 2 )' P= 2- (5) 

Typical values of c are in the range of a few to ten for type I and perhaps a factor of three smaller 
for type II. There is a weak dependence on mass, such that less massive halos have a larger central 
density (e.g., Cole & Lacey 1996; Tormen, Bouchet, k White 1997; Navarro et al. 1996, 1997; Jing 
& Suto 2000). This is understood in general terms as reflecting the mean density at the redshift 
Zf when the halo initially collapsed, 5 ~ (1 + Zf) 3 . For S7 = 1 this is c ~ cr(M), or c oc M - ( 3+n )/ 6 
in a scale-free model. 



2.2. Halo Mass Function 

The number density of halos with mass M within a logarithmic interval is often approximated 
by the prescription of Press & Schechter (1974), 

_dn_ = [2 dlna- 1 ^ 2/2 = ^_ 

cilnM VWlnMM ' cr(M) ' v; 

where <5 C is a parameter characterizing the linear overdensity at the onset of gravitational collapse, 
and a is the linear rms mass fluctuations in spheres of radius R 

* 2 (M) = J™^^P(k)wHkR), (7) 

where W(x) = 3(sin2; — xcosx)/x 3 is the Fourier transform of a real-space tophat window function. 
The mass M is related to R by M = AirpR? /3. For scale free models with a power law initial power 
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spectrum P oc k n , this is a = (M/M„) ( 3 + n )/ 6 . The parameter M* characterizes the mass scale at 
the onset of nonlinearity, <j(M*) = 1, and is related to the nonlinear wavenumber k n \ (defined as 
}J^ l ^itk 2 dkP{k)/{2itf = 1) by 



4npB(n) _ 4irp 3 
3 fc3 3 *' 



(8) 



where 



£(n) 

S (3+n)/3 



(n + 3) / 
Jo 



sin 



. 7T 

(n + 2)^ 



r(n + 2; 



9(2~ n )(3 + n 



3/ (n+3) 



(9) 



(— n)(l — n)(3 — n) 

(defined for —3 < n < 1). Various modifications to the Press-Schechter mass function have been 
suggested (e.g., Sheth & Tormen 1999; Lee & Shandarin 1999; Jenkins et al. 2000) to improve the 
accuracy of the original formula. 



2.3. Halo-Halo Correlations 

Dark matter halos do not cluster in the same way as the mass density field. On large scales, a 
bias parameter b is typically used to quantify this difference. Let £haio(?"; M, M') be the two-point 
correlation function of halos with masses M and M', £iin(f) be the linear correlation function for 
the mass density field, and Phaio and P\\ D be the corresponding power spectra. On large length 
scales, we assume a linear bias and write 

£haio(r; M,M') = 6(M)6(M')6in(r), 

P halo (k;M,M') = b(M)b(M')P lin (k). (10) 

Based on the peak and the Press-Schechter formalism, Mo & White (1996) developed a model for 
the linear bias b(M), which is later modified by Jing (1998) to be 

„2_ n / j n 0.06-0.02 „ . 



6(M » = ( 1+ — JU? +1 J • v= ^ky (11> 

The original formula for b(M) by Mo & White includes only the first factor above; the second 
factor, dependent on the primordial spectral index n, is obtained empirically for an improved fit 
to simulation results at the lower mass end (Jing 1998). In this bias model, b(M) is below unity 
for M < M* (where a{M tf ) = 1) and reaches ~ 0.5 for M < 0.01M*. Small dark matter halos 
are therefore anti-biased relative to the mass density. For M > M*, b{M) increases monotonically 
with the halo mass and reaches b ~ 10 at M ~ 100 M*. Nonlinear effects on the bias have been 
studied (Kravtsov & Klypin 1999 and references therein), but they are unimportant in our model 
because the halo-halo correlation terms contribute significantly only on large length scales in the 
linear regime (see §3 and 4). 
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Similarly, we use higher order bias parameters to relate the higher-order correlation functions 
for halos and mass density. In this paper we examine the three-point correlation function C( r i> r 2 5 r 3) 
and its Fourier transform, the bispectrum B(k±, k 2 , k 3 ) (see §4 for a more detailed discussion). On 
large length scales where the amplitude of 6 is small, perturbation theory can be used to relate the 
lowest order contribution to the bispectrum of the mass density to the linear power spectrum Pi in 
(Fry 1984): 



B (0) (*i, k 2 , h) = F 12 Pi in (fci)P Un (A: 2 ) + F 23 PU(k 2 )PU(k 3 ) + F 31 Pu n {k 3 ) Piin(h) , 

Fij = y + (k i /k j + k j /k i )(k i -k j ) + ^(k i -kjf. (12) 



Using this perturbative result and the results of Mo, Jing, h White (1997), we can write the halo 
bispectrum as 

£haio(&i, &2, k 3 ; M, M', M") = [b(M)b(M')b(M") F 12 + b{M)b(M')b 2 (M")] P hn (h) Py m (k 2 ) 

+ [b(M)b{M')b{M") F 23 + b(M)b 2 {M')b(M")] P hn (k 2 ) P Vm (k 3 ) 
+ [b(M)b{M')b{M") F 3l + b 2 (M)b{M')b{M")] P hn (k 3 ) fi^fa) , 

(13) 

where b(M) is given by equation (11), and the quadratic bias parameter b 2 (M) is 

(u 2 -l) (v\ 2 , 2 



w = 2iV + Ui {v " 3) - (14) 



For the special equilateral case of k\ = k 2 = k 3 = k, equation (13) simplifies to 
BZ (k;M,M',M") = 



y b{M)b(M')b(M") + b{M)b(M')b 2 (M") 



+ b(M)b 2 {M')b(M") + b 2 {M)b(M')b(M")} P£ n (k). (15) 

In practice, the terms involving b 2 (M) in equations (13) and (15) make only a small net contribution. 
For simplicity, we will therefore not include this term in the subsequent derivations and calculations. 



3. Two-Point Statistics: f (r) and P(k) 

We now construct our analytic halo model for the two-point correlation function £(r) and the 
power spectrum P(k). The two-point correlation function of the cosmological mass density field 
5 = Sp/p is 

Z(r) = (6(x)6(x + r)). (16) 

The Fourier transform of £(r) is the mass power spectrum P(k) = J d 3 r e~ tk ' r £(r), which is related 
to the density field in /c-space by (5(ki)6(k 2 )) = P{k\) (27r) 3 <5z)(fci + k 2 ) , where do is the Dirac 
delta-function. 
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The two-point correlation function measures the excess probability above the Poisson distri- 
bution of finding a pair of objects with separation r (Peebles 1980). The objects can be taken to 
be dark matter particles, most of which cluster gravitationally in the form of dark matter halos. 
One should therefore be able to express £ for the density field in terms of properties of dark matter 
halos. In this picture, we can write the contributions to £ as two separate terms, one from particle 
pairs in the same halo, and the other from pairs that reside in two different halos. In realistic 
cosmological models, dark matter halos exhibit a spectrum of masses that can be characterized by 
a distribution function dn/dM, and the halo centers are spatially correlated. Taking these factors 
into consideration, we can write the two-point correlation function for 5 in terms of the halo density 
profile u(x), halo mass function dn/dM, and halo-halo correlation function £halo discussed in §2. 
We write 

£W=£i h (r) + 6y,(r), (17) 

where the subscripts "l/i" and "2/i" denote contributions from particle pairs in "1-halo" and "2- 
halos", respectively, and 

Cifc(r) = J d 3 r' JdM-^5 2 u(r'/R s ) u(\r' + r\/R s ) 

6*(r) = JdZr'dZr"JdM>^S'u(r'/R>)JdM»^5»u^^ 

= Jd 3 r'd 3 r" JdM^6u(r'/R s )b(M) J dM 6 u{r" / R s ) b(M) 

X 6in(|r'-r" + r|). ~ ~ (18) 

These expressions arise from averaging over displacements r', r" of halo centers from the particle 
positions r±, where r = \r\ — r 2 \- In the last expression above, we have used the bias model of 
equation (10) to relate the halo-halo correlation function £halo to the linear correlation function £ii n 
of the mass density field. 

As we will show in §5, the dominant contribution to the two-point correlation function in the 
nonlinear regime on small length scales is from the first, 1-halo term for particle pairs that reside 
in the same halos. This makes intuitive sense, because closely spaced particle pairs are most likely 
to be found in the same halo. This term is determined by the convolution of the dimensionless 
density profile with itself, 

\{x) = j d 3 yu(y)u(\x + y\). (19) 

For many forms of u(x), the angular integration in this equation is analytic, and A can be reduced 
to a simple one-dimensional integral over y. For some special cases, A can even be reduced to an 
analytic expression. We leave the detailed results for A to the Appendix. 

In Ai-space, the convolutions in equation (18) for £(r) become simple products. Using u(q) 
to denote the Fourier transform of u(x), where u(q) = f d 3 xu(x) e~ tq ' x , we can readily transform 
equation (18) into expressions for the mass power spectrum: 

P(k) = P lh {k) + P 2h {k) , (20) 
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where the 1-halo and 2-halo terms are 

?3 x~n.n m2 



Pih(k) = J dM-^[R 3 Ju(kR s 

P2h(k) = y*dM^^^ ( ^ ) y' dM /^ 7 ^3^~ (fcjR / ) p haio(A , ) (21) 



/ dM — 

J dM 



- 2 

R 3 s 6u(kR s )b(M) Pa n (k) . 



To arrive at the last expression above, we have again used the bias model of equation (10). For 
computational efficiency, we find that the algebraic expressions 

47r{ln(e + l/ g )-ln[ln( e +l/g)]/3} 

Ul{(l) = (1+9 1.1)(2/1.1) ' P =1 

47r{ln( e + l/g) + 0.251n[ln(e+l/g)]} 3 
Ull{q) = 1 + 0.8?" ' P= 2 (22) 

provide excellent fits for the profiles of Navarro et al. (1997) and Moore et al. (1999), with less 
than 4% rms error for form I and less than 1% rms error for form II. The functional form is chosen 
to reproduce the asymptotic behaviors: u ~ 47rlng at small q (with no radial cutoff), and u oc q~ 2 
(type I) and u oc q~ 3 / 2 (type II) at large q. 

The two-point £(r) and P(k) can now be computed analytically from equations (18) and (21). 
The inputs are equation (2) or (22) for the halo density profile u(x) or u(q), equations (4) and 
(5) for R s and 5, equation (6) for the halo mass function dn/dM, and equation (10) for the halo- 
halo correlation function. Since the halo density profile appears to have a nearly universal form 
regardless of background cosmology, £(r) and P (k) depend on cosmological parameters mainly 
through <t(M) of equation (7) and the halo concentration c(M) or central density S(M). (See Ma 
& Fry 2000c for a more detailed discussion of c(M).) 



4. Three-Point Statistics: ( and B 

Here we construct our analytic halo model for the three-point correlation function £ and the 
bispectrum B. The joint probability of finding three objects in volume elements dV±,dV2, and dV^ 
is given by 

dP = [1 + £( ri ) + £(r 2 ) + £(r 3 ) + ((n,r 2 ,r 3 )} n 3 dV 1 dV 2 dV 3 , (23) 

where £(r) and C( r i> r 2; r 3) are the two- and three-point correlation functions, respectively, n is 
the mean number density of objects, and r±,r 2 and r 3 are the lengths of the sides of the triangle 
defined by the three objects (Peebles 1980). The Fourier transform of the three-point correlation 
function C( r i 5 r 2, ^3) is the bispectrum B(k\,k 2 , k 3 ), which is related to the density field in /c-space 
by ( S(ki)S(k 2 )S(ks) ) = B(k\, k 2 , ks) (2ir) 3 5D(ki + k 2 + fc 3 ) . The bispectrum depends on any three 
parameters that define a triangle in /c-space. A particular simple configuration to study is the 
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equilateral triangle {k\ = k 2 = k 3 = k), and in this case the bispectrum B eq depends only on a 
single wavenumber. 

Similar to the two-point halo model of §3, we can write the contributions to the three-point 
correlation function £ of the mass density as three separate terms, each term representing particle 
triplets that reside in a single halo, two distinct halos, or three distinct halos. Taking into account 
the halo mass distribution and halo-halo correlations discussed in §2, we obtain 

((ri,r 2 ,r 3 ) = (i h (ri,r 2 ,r 3 ) + ( 2h (n, r 2 , r 3 ) + (3h(n,r 2 ,r 3 ) , (24) 

where the separate 1-halo, 2-halo, and 3-halo terms are 

/f dn - 

d 3 r / dM — 5 3 u(r/R s ) u(\r + n - r 2 \/R s ) u(\r + n - r 3 \/R s ) 

C2ft(ri,r 2 ,r 3 ) = J d\ d\> J dM 5 2 u(r/R s ) u(\r + n - r 2 \/R s ) J dM> 5> u(r> / R> s ) 

x Chaio(k - r + n - r 3 |) + sym.(l,2,3) (25) 
C3fc(ri,r 2 ,r 3 ) = J d 3 rd 3 r'd 3 r" J dM^8u(r/R s ) J dM' ^ 6> u(r' / R' s ) 



J dM" 



flU — 

' 5"u{r"/R'^) ( halo (r + r 1 ,r' + r 2 ,r" + r 3 ) 



dM" 



The dominant contribution to the three-point correlation function in the nonlinear regime is from 
the first term which comes from particle triplets that reside in the same halo. This term is 
determined by the convolution 7(0:1, x 2 ) = f d 3 yu(y) u(\y + x\\) u(\y + x 2 \) of three factors of the 
density profile u(x), and is analogous to the convolution A in equation (19) for the one-halo term 
£,ih in the two-point correlation function. 

The bispectrum of the mass density field 5 in £>space can be obtained by Fourier transforming 
the equations above. We find 

B(ki, k 2 , k 3 ) = B lh {ki, k 2 , k 3 ) + B 2h (ki, k 2 , k 3 ) + B 3h {ki, k 2 , k 3 ) , (26) 



where 



B lh (k 1 ,k 2 ,k 3 ) = J dM ^ [ r3 s $ u(kiR s )] [R 3 6 u(k 2 R s )] [R 3 S 5 u(k 3 R s )] 
B 2h (h,k 2 ,k 3 ) = JdM^i&JuikMU&JufaRs)] 

x f dM '-^p R f & u(k 3 R' s ) P halo (fc 3 ; M, M') + sym.(l,2,3) (27) 
B 3h (h,k 2 ,k 3 ) = j 'dM^RlduihRs) J ' d M'-^R' 3 S'u(k 2 R' s ) 

x I dM" -^Rfd"u(k 3 R'l)B halo (k 1 ,k 2 ,k 3 ;M,M',M"). 
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The halo-halo power spectrum Phaio(^) an d bispectrum -Bhaio(^i) ^2, ^3) are related to the linear 
mass power spectrum Pn n (k) by equations (10) and (13). 

The expressions for the mass bispectrum above simplify considerably for the equilateral triangle 



configuration, and 
where 



(28) 



BH(k) = JdM^[Rl5u(kR 
KW = 3 



B2(k) = 



dM^[Rl5u(kR s )fb{M) 
dn 



j "dM-^R 3 s 5 u{kR s ) b(M) P lin (k) (29) 



dM ^R 3 Ju(kR s )b(M) 



3 12 



Here we have written out explicitly the bias factors b(M) using equations (10) and (15), and we 
have neglected terms with 62 (M) as discussed in §2.3. 



5. iV-body Experiments and Numerical Results 

In this section we compare the predictions of our analytical model described in §2, 3, and 4 with 
results from cosmological iV-body simulations. We examine two cosmological models: an n = — 2 
scale-free model and a low-density ACDM model. These are the same simulations studied in Ma 
& Fry (2000a). The n = —1 simulation has 256 3 particles and a Plummer force softening length of 
L/5120, where L is the box length. The ACDM model is spatially flat with matter density £l m = 0.3 
and cosmological constant S7a = 0.7. This run has 128 3 particles and is performed in a (100 Mpc) 3 
comoving box with a comoving force softening length of 50kpc for Hubble parameter h = 0.75. 
The baryon fraction is set to zero for simplicity. The primordial power spectrum has a spectral 
index of n = 1, and the density fluctuations are drawn from a random Gaussian distribution. 
The gravitational forces are computed with a particle-particle particle-mesh (P 3 M) code (Ferrell 
&; Bertschinger 1994). We compute the density field 5 on a grid from particle positions using 
the second-order triangular-shaped cloud (TSC) interpolation scheme. A fast Fourier transform is 
then used to obtain 5 in /c-space. The /c-space TSC window function is deconvolved to correct for 
smearing in real space due to the interpolation, and shot noise terms are subtracted to correct for 
discreteness effects. We then compute the second and third moments of the density amplitudes in 
Fourier space. 

We show results for the power spectrum as the dimensionless variance A(k) = 4-7r/c 3 P(/c) /(2-7r) 3 . 
A useful dimensionless three-point statistic is the hierarchical three-point amplitude 

nl , , , x B(k l ,k 2 ,k 3 ) 

to, M - p{ki)p{k2) + p {k2)P{k3) + P^Pih) ■ ^ 
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The three-point amplitude Q has the convenient feature that for the lowest nonvanishing result in 
perturbation theory, Q is independent of time and the overall amplitude of P; for scale-free models 
with a power-law P, Q is independent of overall scale as well. To lowest order, it follows from 
equation (12) that the equilateral bispectrum has a particularly simple form, B(°\k) = ^ ^u n (^) > 
and we have (k) = | , independent of the power spectrum. 

5.1. Synthetic Halo Replacement 

To investigate the numerical effects of limited resolution in the simulations, we have experi- 
mented with the distribution of matter in halos identified in the simulations. In these experiments, 
we keep the locations and masses of the halos unchanged but redistribute the subset of particles 
which lies within the virial radius -R200 (the radius within which the mean overdensity is 200) of 
each halo according to a prescribed density profile. We then recompute the two- and three-point 
statistics A and Q from the redistributed particle positions as well as the original non-halo parti- 
cles, which remain at their original positions. By using density profiles obtained empirically from 
higher-resolution simulations of individual halos, this recipe allows us to model accurately the inner 
regions of the halos on scales below the numerical softening length scale while at the same time 
preserving all the large-scale information available in the large parent simulation. This technique 
should also be useful for other studies that are sensitive to the inner halo density profiles, for 
example the ray-tracing method in gravitational lensing. 

Ma & Fry (2000a) have used this replacement technique to experiment with synthetic halos 
that follow a pure power-law profile u oc r~ e . It is found that A(k) and Q(k) at high-/c indeed 
obey A(k) oc k 2e ~ 3 and Q(k) oc k 3 ~ e as predicted by the simple power-law model of Peebles (1974). 
The scaling works even in the presence of the full distribution of matter outside the halo cores. 
Here we extend this replacement technique to more realistic halo profiles of equation (2). Figures 
1 and 2 illustrate the effects on the matter power spectrum and bispectrum when the original 
halos in large cosmological simulations are replaced by synthetic halos with the density profile 
uu = l/(x 3 / 2 + x 3 ) of equation (2). For the n = —2 scale-free model, the concentration parameter 
is taken to be c(M) = 3(M^/M) 1 / 6 , which is consistent with Navarro et al. (1997) and has the 
expected scaling with mass, c oc M~( 3+n )/ 6 , in a scale-free model. For ACDM models, we use 
c(M) = 5(M*/M) 1 / 6 as suggested by Figure 3 of Moore et al. (1999). We note, however, that 
c(M) from various recent simulations has shown a large scatter, and its functional form depends 
on the exact form of the density profile used. For the ACDM model and form ujj = l/(x 3 / 2 + x 3 ), 
for example, a flatter and smaller c(M) = 3(M ! „/M) a084 appears to be preferred by Jing & Suto 
(2000) and Navarro et al. (1997). The results of Tormen et al. (1997) and Cole & Lacey (1996) 
are also only marginally consistent with each other. A more detailed investigation of the different 
forms of c(M) can be found in Ma k Fry (2000c). 

In Figures 1 and 2, the agreement at low values of k between the original and synthetic halos 
is excellent, confirming that the correlation functions on larger length scales are insensitive to 



-12- 



the spatial distribution of particles in the halo cores. The only significant difference between the 
simulation and synthetic halos appears at small length scales, where the coarser resolution of the 
simulation blurs out the structure of the inner halo and results in an inner profile flatter than in 
equation (2). This effect is manifested in the bending over of the dashed curves for P(k) in Figures 
1 and 2 at high k, and is corrected for when the synthetic halos are used. 

5.2. TV-body Results vs. Analytic Halo Model 

We now proceed to compare the predictions of the analytic model of §2 - §4 with the numerical 
results from cosmological simulations. Figures 3 and 4 show the fc-space density variance A(k) 
(upper panel) and the three-point amplitude Q C q(k) for equilateral triangles for the n = —2 scale- 
free model and the ACDM model. The solid black curves are the model predictions computed 
from equations (21) and (29). The contribution from the single-halo and multiple-halo terms are 
shown separately as dashed curves. For the density profile, we use the same uji = l/(x 3 / 2 + x 3 ) 
and concentration parameters as in Figures 1 and 2. For the mass function, we use the Press- 
Schechter formula but reduce its overall amplitude by 25%, which we find necessary in order to 
match the halo mass functions for our numerical simulations. This overestimation of halo numbers 
with M ~ M* by Press-Schechter is a well known result reported in many other studies (see Jenkins 
et al. 2000 and references therein). The mass limits for the integrals in equations (21) and (29) 
do not significantly affect the model predictions for the total A or Q. Raising the lower mass limit 
does reduce the contribution from lower mass halos and hence lower the high-fc amplitudes of the 
multiple halo terms A 2 h, Q2h, an d Q3h, but these terms make negligible contributions to the total 
A and Q. 

As discussed in §3 and 4, the nonlinear parts of both the two- and three-point statistics are 
determined by the dominant 1-halo term because the closely spaced particle pairs and triplets 
mostly reside in the same halos. The multiple-halo terms are therefore significant only on larger 
length scales comparable to the separation between halos. Their inclusion, however, is necessary 
for the transition into the linear regime. 

For the n = —2 model in Figure 3, we plot the results against the scaled k/k n \, where k n \ 
characterizes the length scale that is becoming nonlinear and is defined by f Q nl d 3 k P\m(a, k) / (2-7r) 3 = 
1. Three time outputs are shown, where the expansion factor (1 initially) and k n \ (in units of 
2vr/L) are: (a, fc n i) = (13.45, 29), (19.03, 14.5), and (26.91,7.25) (from left to right). For the two- 
point A(k), the agreement between the model prediction and the simulations is excellent. The 
three simulation outputs also overlap well, indicating that self-similarity is obeyed, as reported in 
Jain & Bertschinger (1998). For the three-point Q e q, however, self-similar scaling does not hold 
as rigorously (Ma & Fry 2000a). It is interesting to note that the analytic prediction agrees most 
closely with the earliest output (a,k n i) = (13.45,29) (green curve). This provides further evidence 
to the suggestion of Ma Sz Fry (2000a) that the later outputs of the n = —2 simulation may 
be affected by the finite volume of the simulation box. For the ACDM model in Figure 4, the 
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analytic model again provides a good match to the iV-body results within the fluctuations among 
the simulations. We illustrate the numerical effects due to box sizes by showing results from two 
runs with volume (100 Mpc) 3 and (640 Mpc) 3 . The model predictions extend well beyond the 
resolution of the simulations. 

The real-space two-point correlation function for the n = —2 and ACDM models is shown in 
Figures 5 and 6. For the halo model predictions, we have chosen to show only the results for the 

1- halo term ^ because this term dominates the interesting nonlinear portion of £. The agreement 
between the halo model (dashed curves) and the simulations (symbols) is again excellent. For the 

2- halo terms ^2hi the computation can be done more easily in fc-space as shown in Figures 3 and 4, 
so we do not include them here. 

For comparison, we plot in Figures 3-6 the results from the commonly used fitting formulas for 
the nonlinear power spectrum (Hamilton et al. 1991; Jain et al. 1995; Peacock & Dodds 1996; Ma 
1998; Ma et al. 1999). While the formulas provide a good approximation to A(k) up to k/k n \ ~ 50 
for the n = —2 model and k ~ 20 h Mpc -1 for the ACDM model, the figures show that significant 
deviations occur at higher k, and the fitting formula and our current model predict different high-fc 
slopes for A(k). Since the high-A; behavior of the fitting formulas has been constructed to obey the 
stable clustering prediction, this discrepancy has an important implication for the validity of stable 
clustering, which we discuss briefly in the next section and at length in Ma & Fry (2000c). 

6. Discussion 

We have constructed a physical model for the correlation functions of the mass density field 
in which the correlations are derived from properties of dark matter halos. We have described in 
detail the input, construction, and results of this model in §2 - §5. We now examine more closely 
its physical meanings and implications in three separate regimes. 

On scales larger than the size of the largest halo, the contributions from separate halos dom- 
inate, and (by design) the model reproduces the results of perturbation theory. On intermediate 
scales, 1/R* < k < 1/R S (M*), because of the exponential cutoff in the mass function dn/dM 
at the high mass end, the contribution to the volume integrals in equation (18) is dominated by 
the large-r regime where the halo profiles are roughly r -3 . The correlation functions therefore 
behave approximately as predicted by the power-law model with e = 3, i.e., A oc k 2 "- 3 ~ k s and 
Q oc k 3 ~ e ~ constant. This is why Q exhibits an approximately flat plateau at intermediate k in 
the bottom panels of Figures 3 and 4. 

On the smallest and most nonlinear scales, the correlation functions probe the innermost 
regions of the halos. Intriguingly, the halo model predicts on these scales a behavior that is different 
from either the frequently-assumed stable clustering result of A(k) oc /c 7 with 7 = (9 + 3n) /(5 + n) 
(Davis & Peebles 1977), or the power-law profile result of 7 = 2e — 3. The implication of departure 
from stable clustering is significant because all the fitting formulas for the nonlinear P(k) in the 
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literature (see §5.2) have been constructed to approach the stable clustering limit at high k. A 
more detailed study on the criteria for stable clustering in this model is given in a separate paper 
(Ma & Fry 2000c). 

The origin of the deviation from stable clustering in the model at high-/c can be understood as 
follows. For the two-point function, as k becomes large, the one-halo integral P\h{k) in equation (21) 
converges before the exponential cutoff, and is dominated by contributions near the mass scale for 
which kR s = 1. The behavior now depends on the mass distribution function. The various mass 
functions discussed in §2.2 have the same general behavior of dn/dM oc M~ 2 v a e~ u2 / 2 , where 
v = 5 c /a. The Press-Schechter form assumes a = 1 (see eq. [6]), while others (e.g., Sheth & 
Tormen 1999; Jenkins et al. 2000) suggest a flatter slope of a ~ 0.4 for the lower mass halos. Since 
the scale radius R s depends on mass as R s = R 2 oo/c oc M 1 / 3 /M _ ( 3+n )/ 6 oc M^ 11 ^ 6 , and R 3 S 8 oc M 
(up to logarithmic factors), we find from equation (21) that the power spectrum at high k goes as 

A(k) « Ai h (ife) oc fc 3 J dM v a u 2 (kR s ) . (31) 
Changing variables to y = kR s oc k (M/M*)( 5+n )/ 6 , we see that 

7-(t£)-(§s). 

where the first term in 7 is the prediction of stable clustering. The departure arises from the factor 
v a in the mass function, and would vanish only if a = or n = —3. This is the origin of the 
difference in A(fc) at high k between the model prediction (solid curves) and the fitting formula 
(dotted curves) shown in Figures 3 and 4. 

For the three-point function, the one-halo integral B^ in equation (29) converges (barely, for 
p = I and n = —2), giving 

*w K 73 ^(^)- Q (|^) 

Q cq (k) oc fc a ( 3+n ^ 5+n ) (33) 

This again disagrees with the prediction of stable clustering that Q is constant, but it appears to 
be consistent with numerical simulations as shown in Figures 3 and 4. 

For yet higher order correlations, details of the halo profile begin to matter. For p = 1, the 
pattern of equations (32) and (33) persists to all orders, but for p = | they apply only for the 
two- and three-point functions; for four-point and higher functions the nonlinear scale M* and 
7 n = n p — 3 for n > 4. Thus there seems to be some potentially interesting behavior that is tested 
only in the four-point function and higher. 



(32) 
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7. Summary 

We have presented an analytic model for the two- and three-point correlation functions £(r) 
and C( r i, f2, rs) of the cosmological mass density field and their Fourier transforms, the mass power 
spectrum P(k) and the bispectrum B(ki,k2,ks). In this model, the clustering statistics of the 
density field are derived from a superposition of dark matter halos with a given set of input halo 
properties. These input ingredients include realistic halo density profiles of equation (2), halo mass 
distribution of equation (6), and halo-halo spatial correlations of equations (10) and (15). The 
main results of the model are given by equations (18) and (21) for the two-point statistics £ and 
P, and by equations (25) and (29) for the three-point statistics £ and B. This model provides a 
rapid way to compute the correlation functions over all length scales where the model inputs are 
valid; it also gives a physical interpretation of the clustering process of matter in the universe. 

We have tested the validity of this model by comparing its predictions with results from 
cosmological simulations of an n = —2 scale-free model and a ACDM model. As Figures 3-6 
illustrate, the model describes well the simulation results spanning the entire range of behavior 
from the perturbative regime on large scales to the strongly nonlinear regime on small scales. To 
probe the critical high-/c range in the deeply nonlinear regime, we have used a halo replacement 
technique to increase the resolution of the large parent simulations. As Figures 1 and 2 illustrate, 
this method of replacing the original halos that suffer from numerically softened cores with synthetic 
halos of analytic profiles is a reasonable way to improve the resolution of numerical simulations. By 
using density profiles obtained empirically from higher-resolution simulations of individual halos, 
this recipe allows us to model accurately the inner regions of the halos on scales below the numerical 
softening length scale, while at the same time preserving all the large-scale information available 
in the large parent simulation. This technique should also be useful for other studies that depend 
on the inner halo density profiles, for example, the ray-tracing method in gravitational lensing. 

Given that dark matter halos in simulations (and presumably in nature) are not perfectly 
spherical, cleanly delineated objects, it is intriguing that the model constructed in this paper works 
as well as it does at matching the simulation results. Nevertheless, this analytic model provides 
a good qualitative and quantitative description over the entire range of scales covered by the 
simulation, and it can be used to make predictions beyond these scales. This is the first model 
prescription that successfully reproduces both two- and three-point mass correlations. We believe 
that it will prove to be a generally useful framework. 

We have enjoyed stimulating discussions with John Peacock and David Weinberg. We thank 
Edmund Bertschinger for valuable comments and for providing the n = — 2 scale-free simulation. 
Computing time for this work is provided by the National Scalable Cluster Project and the Intel 
Eniac2000 Project at the University of Pennsylvania. C.-P. M. acknowledges support of an Alfred 
P. Sloan Foundation Fellowship, a Cottrell Scholars Award from the Research Corporation, a Penn 
Research Foundation Award, and NSF grant AST 9973461. 
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A. Appendix 



In this Appendix we display analytic forms for the convolution of the dimensionless profile 
shape 

A(x) = [ d 3 yu(y)u(\x + y\) (Al) 



discussed in §3. These analytic expressions are useful for computing the nonlinear two-point correla- 
tion function £ of the mass density field, which is dominated by the 1-halo term in equation (18) 
and is related to A by 



(r) = / 



dM^5 2 ^ A(r/jRs)) for£>l. 



(A2) 



For the type-I profile uj of equation (2), the angular integration in equation (Al) is analytic, 
and A is reduced to a simple integral 



A/(x) 



2tt 



ydy 



(2-p)xJ yP{l + yf-P 



(x + y) 



2-p 



\x-y 



2-p 



{l + x + yf-P (1 + \x - y\) 2 - p _ 



For the special case p = 1, this integral can be further reduced to the analytical form 



A/(x) = 



8tt 



(x 2 + 2x + 2) ln(l + x) 
x (x + 2) _ 



x 2 (x + 2) 

For n// of equation (2), we are able to simplify A to 

2tt f 00 ydy 



- 1 



p=l. 



\ H (x) 



x Jo y p (i + y 3 ~ p ) 



F P (x,y) , 



(A3) 



(A4) 



(A5) 



where the function F p (x,y) represents the angular part of the integration in equation (Al) and 



F p {x,y) 



x+y 



z dz 



x—y\ 



ZP(1 + Z 3 ~P) 



(A6) 



The integral in F p can be reduced to analytic forms for special values of p. Here we display the six 
cases p = 0, 1/2, 1, 3/2, 2, and 5/2: 



- |2\/3tan _1 
6 



-l + 2(g + y) 
V3 



+ ln 



1 — (x + y) + (x + y) 2 
l + 2(x + y) + (x + y) 2 _ 



'1/2 



- {replace (x + y) above with |x — y|} 
1 I -2,/^ tan- ( 1 + ^- 4 ^ 



(A7) 
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F, 



3/2 



"5/2 



+41n(l + V^Ty) -(1 + V5)ln 



! + -(-! + >/5)VaT+y + x + y 



-(1 - >/5) In 



1 



— — {replace (x + y) above with \x — y\} 



tan 1 (x + y) — tan 1 (\x — y\) 
1 



2^ tan 



-i 



In 



x + y 
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2 
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\x-y\ 



+ In 
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y/\x-y\ Vx + y 
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Fig. 1. — Effects on the power spectrum (upper panel) and bispectrum (lower panel) when the dark 
matter halos in an n = -2 scale- free simulation are replaced with synthetic halos of density profile 
uii{x) = l/(x 3 / 2 + x 3 ) and concentration parameter c(M) = 3(M*/M) 1 / 6 (see §2 for definitions). 
The dashed and solid curves are for the original and the redistributed particles, respectively. They 
agree up to k ps 200 (in 2ir/L), beyond which the dashed curves deviate due to the finite resolution 
in the original simulation. The dotted curves show the linear A and the nonlinear fitting of Jain et 
al. (1995) in the upper panel, and the lowest-order perturbative result = 4/7 in the bottom 
panel. 
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Fig. 2. — Same as Fig. 1 but for a low-density CDM simulation with 0, m = 0.3, f^A = 0.7. The 
synthetic halos have the un(x) = l/(x 3 / 2 + x 3 ) profile and concentration parameter c(M) = 
5(M*/M) 1 /6. A gain, the original (dashed) and redistributed (solid) particles have similar A(k) 
and Qeq(k) up to the simulation resolution of k « 20 h Mpc -1 . 



-22- 




Fig. 3. — iV-body results vs. predictions of the analytic model of §2-§4 for the power spectrum 
(upper) and bispectrum (lower) for the n = —2 scale-free model. The dashed curves show the 
separate contributions to A and Q e q computed from the single- and multiple-halo terms of eqs. (21) 
and (29); the solid black curves show the sum predicted by the model. The colored curves show 
the TV-body results, where synthetic halos have been used to extend the curves to higher k. (The 
same density profile and c(M) are used for the synthetic halos and the analytic model.) Three 
simulation outputs are shown, where the expansion factor (1 initially) and nonlinear wavenumber 
(in units of 2vr/L) are: (a,k nl ) = (13.45, 29), (19.03, 14.5), and (26.91, 7.25) (from left to right in 
green, blue and red). Three of the four dotted curves are the same as in Fig. 1; the rising one in 
the bottom panel shows the 1-loop Q. 
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Fig. 4. — Same as Fig. 3 but for the low-density CDM simulation with f2 m = 0.3, = 0.7. The 
red and green curves are computed from a (100 Mpc) 3 and a (640 Mpc) 3 simulation, respectively. 
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Fig. 5. — iV-body results vs. predictions of the analytic model of §2-4 for the two-point corre- 
lation function £(r) for the n = —2 model. The dashed curve shows the 1-halo term £ih(r) of 
equation (18) from our analytic model. The solid squares show £(r) computed directly from an 
iV-body simulation. The two agree very well for r/r Q \ < 1. The dotted curve shows the linear 
theory £i in (r) = r n i/r. 
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Fig. 6. — Same as Fig. 5 but for the ACDM model. The symbols show £(r) computed from a (100 
Mpc) 3 (open circles) and a (640 Mpc) 3 (solid squares) iV-body simulation. The dotted curves show 
£iin(f) from the linear theory (lower curve) and the nonlinear £(r) (upper curve) given by the fitting 
formula of Ma (1998). 



